home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / sys / invhyp.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-03-15  |  2.3 KB  |  109 lines
  1. /* sys/invhyp.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_machine.h>
  24.  
  25. double gsl_acosh (const double x);
  26.  
  27. double
  28. gsl_acosh (const double x)
  29. {
  30.   if (x > 1.0 / GSL_SQRT_DBL_EPSILON)
  31.     {
  32.       return log (x) + M_LN2;
  33.     }
  34.   else if (x > 2)
  35.     {
  36.       return log (2 * x - 1 / (sqrt (x * x - 1) + x));
  37.     }
  38.   else if (x > 1)
  39.     {
  40.       double t = x - 1;
  41.       return log1p (t + sqrt (2 * t + t * t));
  42.     }
  43.   else if (x == 1)
  44.     {
  45.       return 0;
  46.     }
  47.   else
  48.     {
  49.       return GSL_NAN;
  50.     }
  51. }
  52.  
  53. double gsl_asinh (const double x);
  54.  
  55. double
  56. gsl_asinh (const double x)
  57. {
  58.   double a = fabs (x);
  59.   double s = (x < 0) ? -1 : 1;
  60.  
  61.   if (a > 1 / GSL_SQRT_DBL_EPSILON)
  62.     {
  63.       return s * (log (a) + M_LN2);
  64.     }
  65.   else if (a > 2)
  66.     {
  67.       return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
  68.     }
  69.   else if (a > GSL_SQRT_DBL_EPSILON)
  70.     {
  71.       double a2 = a * a;
  72.       return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
  73.     }
  74.   else
  75.     {
  76.       return x;
  77.     }
  78. }
  79.  
  80. double gsl_atanh (const double x);
  81.  
  82. double
  83. gsl_atanh (const double x)
  84. {
  85.   double a = fabs (x);
  86.   double s = (x < 0) ? -1 : 1;
  87.  
  88.   if (a > 1)
  89.     {
  90.       return GSL_NAN;
  91.     }
  92.   else if (a == 1)
  93.     {
  94.       return (x < 0) ? GSL_NEGINF : GSL_POSINF;
  95.     }
  96.   else if (a >= 0.5)
  97.     {
  98.       return s * 0.5 * log1p (2 * a / (1 - a));
  99.     }
  100.   else if (a > GSL_DBL_EPSILON)
  101.     {
  102.       return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
  103.     }
  104.   else
  105.     {
  106.       return x;
  107.     }
  108. }
  109.